/***
   *
   * This program uses the discrete cosine transform to compress a pixel 
   * image by a factor of 4:1 while preserving its information content.
   * Block encoding is used with a block size of 8 x 8.
   *
   * This program is based on the routines and algorithms found in the book
   * "C Language Algorithms for Digital Signal Processing" By P. M. Embree
   * and B. Kimble.
   *
   * Copyright (c) 1992 -- Mazen A.R. Saghir -- University of Toronto   
   * Modified to use arrays - SMP 
   * Modified by Catherine Wong, 1998
   *    - deals with .pgm image files now, instead of text files
   *    - combined with expand.c
   *
 ***/


#ifndef __COMPRESS_H
#define __COMPRESS_H

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>


#define  ROUND(f)          (((f)<0.0) ? (int)((f)-0.5) : (int)((f)+0.5))
#define  B                 8


void image_dct(unsigned char *in, unsigned char *out, int inrow, int incol);
				/* Image Discrete Cosine Transform */ 

void block_dct(float *block, float cos1[B][B], float cos2[B][B]);
				/* Block Discrete Cosine Transform */ 

unsigned char pack2in8(int a, int b);   /* Packs 2 ints into 1 unsigned char */
void init_cos_matrices (float cos1[B][B], float cos2[B][B]);

#endif

void init_cos_matrices (float cos1[B][B], float cos2[B][B]);

main( int argc, char** argv )
{
   unsigned char *in;           /* Points to input image */ 
   unsigned char *out;          /* Points to output compressed */ 
   FILE *infp, *outfp;          /* Original and compressed image file ptrs */ 

   int inrow, incol;


  /*
   * Check argument count and process arguments 
   */

/*   if (argc != 3) { 
      fprintf(stderr, "Usage: Compress <in.pgm> <out.z>\n");
      return -1;
   }*/
/*   infp  = fopen(argv[1], "rb");*/
   infp  = fopen("bird.pgm", "rb");
/*   outfp = fopen(argv[2], "wb"); */
   outfp = fopen("comp.out", "wb"); 
   if (infp == NULL || outfp == NULL) return -1;

  /*
   * Read the image 
   */

   if (fscanf(infp,"P5\n%*s%*s%*s%*s%*s%*s%*s\n%d%d\n255\n",&incol,&inrow) != 2)
      return -1;

   in     = malloc(sizeof(unsigned char) * incol * inrow);
   out    = malloc(sizeof(unsigned char) * inrow * incol);
   assert(in != NULL & out != NULL);

   assert(fread(in, sizeof(unsigned char), incol * inrow, infp));
   fclose(infp);

  /*
   * Apply DCT compression to the input image
   */

   image_dct(in,out,inrow,incol);

  /*
   * Write out the compressed image 
   */

   fprintf(outfp, "%d %d\n", 2 * incol, inrow / 8); 
   fwrite (out, sizeof(unsigned char), incol * inrow / 4, outfp);
   fprintf(outfp, "\n");
   fclose(outfp);

   exit(0);
}



/***
   *
   * This function packs 2 4B integers into 1 1B unsigned char
   *
 ***/
 
unsigned char pack2in8(int a, int b)
{
   unsigned char c;
   
   if( a >  7 ) a =  7;
   if( b >  7 ) b =  7;
   if( a < -8 ) a = -8;
   if( b < -8 ) b = -8;

   a &= 0x0f;
   b &= 0x0f;

   c = (unsigned char)( (a << 4) | b );
   return(c);
}


/*
 * Initialize the cosine matrices. "cos2" is the transpose of "cos1" 
 */

void init_cos_matrices( float cos1[B][B], float cos2[B][B] )
{
   int i, j;
   float temp_cos;
   float factor1, factor2;

   factor1 = 2.0 * atan(1.0)/B;
   factor2 = 0.0;
   for (i = 0; i < B; i++) {
      for (j = 0; j < B; j++) {
         temp_cos = cos(factor2 * (2*j + 1))/B;
         cos1[i][j] = temp_cos;
         cos2[j][i] = temp_cos;
      }
      factor2 += factor1;
   }
}
/***
   *
   * This file contains the vectorizable functions
   *
 ***/


  /*
   * Each image is assumed to have a resolution of inrow x incol, with 256 gray
   * levels. Since each block is 8 x 8, (inrow x incol/64) iterations are 
   * needed to compress the image. This is the main loop. 
   * 
   * Images need to be stored as signed characters so that -1 does not turn
   * into 255. 
   */
 
void image_dct(unsigned char *in, unsigned char *out, int inrow, int incol)
{
   int row;       /* General variable, used to simplify row manipulations */ 
   int blkrow;    /* Number of rows of blocks in the input image */
   int vblk;      /* Vertical number of the current block being accessed */

   int col;       /* General variable, used to simplify column manipulations */ 
   int blkcol;    /* Number of columns of blocks in the input image */
   int hblk;      /* Horizontal number of the current block being accessed */

   int i, j;      /* General counters */

   float cos1[B][B];	/* lookup table of cosine values */
   float cos2[B][B];	/* transpose of cos1 array */
   float *block;
   unsigned char *outptr;  /* Used to simplify manipulations of output array */


   init_cos_matrices(cos1,cos2);

   blkrow = inrow/8;
   blkcol = incol/8;

   block  = malloc(sizeof(float) * B * B);
   assert(block != NULL);

   for (vblk = 0; vblk < blkrow; vblk++) {
      for(hblk = 0; hblk < blkcol; hblk++) {

        /*
         * Read next image block. 
         */

         for (i = 0; i < B; i++) {
            for (j = 0; j < B; j++) {
               row = B*vblk + i; 
               col = B*hblk + j;
               (block+B*i)[j] = (float)(in + row*incol)[col];
            }
         }


        /*
         * Find its DCT 
         */

         block_dct(block, cos1, cos2);


        /*
         * Select coefficients, scale, and pack 
         */

/*#define RNDBLK(i,j,p)	(ROUND((block+(B*i+j)*P)[p]))*/

#define BLK(i,j)		((block+B*i)[j])
#define RNDBLK(i,j)		(ROUND(BLK(i,j)))
#define RNDBLKDIV(i,j,k)	(ROUND((BLK(i,j))/k))

         outptr = out + 16*(vblk*blkcol + hblk);
         outptr[0] = (unsigned char)RNDBLK(0,0);
         outptr[1] = (signed char)RNDBLK(0,1);
	 outptr[2] = pack2in8(RNDBLKDIV(0,2,3),RNDBLKDIV(0,3,3));
	 outptr[3] = pack2in8(RNDBLKDIV(0,4,2),RNDBLKDIV(0,5,2)); 
         outptr[4] = pack2in8(RNDBLK   (0,6),  RNDBLK   (0,7)); 
         outptr[5] = (signed char)RNDBLK(1,0);
         outptr[6] = (signed char)RNDBLK(1,1); 
	 outptr[7] = pack2in8(RNDBLKDIV(1,2,2),RNDBLK   (1,3)); 
         outptr[8] = pack2in8(RNDBLK   (1,4),  RNDBLK   (1,5)); 
	 outptr[9] = pack2in8(RNDBLKDIV(2,0,3),RNDBLKDIV(2,1,2));
	 outptr[10]= pack2in8(RNDBLKDIV(2,2,2),RNDBLK   (2,3)); 
	 outptr[11]= pack2in8(RNDBLKDIV(3,0,3),RNDBLK   (3,1)); 
         outptr[12]= pack2in8(RNDBLK   (3,2),  RNDBLK   (3,3)); 
	 outptr[13]= pack2in8(RNDBLKDIV(4,0,2),RNDBLK   (4,1)); 
	 outptr[14]= pack2in8(RNDBLKDIV(5,0,2),RNDBLK   (5,1)); 
	 outptr[15]= pack2in8(RNDBLK   (6,0),  RNDBLK   (7,0)); 
      }
   }

}


/***
  * This function calculates the discrete cosine transform of an image block.
  * The out is returned in the same matrix as the original block. 
 ***/

void block_dct( float *block, float cos1[B][B], float cos2[B][B] )
{
   int i, j, k;
   float sum;

   float tmpblk[B][B];

  /*
   * Multiply the input image block with the cos2 matrix; store the out
   * in the temporary matrix "tmpblk". 
   */

   for (i = 0; i < B; i++) {
      for (j = 0; j < B; j++) {
         sum = 0.0;
         for (k = 0; k < B; k++) {
            sum += (block + B*i)[k] * cos2[k][j];
         }
         tmpblk[i][j] = sum;
      }
   }


  /*
   * Multiply the cosine matrix by the temporary matrix; store the
   * out back in the original matrix.   
   */

   for (i = 0; i < B; i++) {
      for (j = 0; j < B; j++) {
         sum = 0.0;
         for (k = 0; k < B; k++) {
            sum += cos1[i][k] * tmpblk[k][j] ;
         }
         (block + B*i)[j] = sum;
      }
   }
}


